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ABSTRACT 

The gamma-ray burst (GRB) rate is essential for revealing the connection 
between GRBs, supernovae and stellar evolution. Additionally, the GRB rate 
at high redshift provides a strong probe of star formation history in the early 
universe. While hundreds of GRBs are observed by Swift, it remains difficult to 
determine the intrinsic GRB rate due to the complex trigger algorithm of Swift. 
Current studies of the GRB rate usually approximate the Swift trigger algorithm 
by a single detection threshold. However, unlike the previously flown GRB in- 
struments, Swift has over 500 trigger criteria based on photon count rate and 
additional image threshold for localization. To investigate possible systematic 
biases and explore the intrinsic GRB properties, we develop a program that is 
capable of simulating all the rate trigger criteria and mimicking the image thresh- 
old. Our simulations show that adopting the complex trigger algorithm of Swift 
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increases the detection rate of dim bursts. As a result, our simulations suggest 
bursts need to be dimmer than previously expected to avoid over-producing the 
number of detections and to match with Swift observations. Moreover, our re- 
sults indicate that these dim bursts are more likely to be high redshift events 
than low- luminosity GRBs. This would imply an even higher cosmic GRB rate 
at large redshifts than previous expectations based on star-formation rate mea- 
surements, unless other factors, such as the luminosity evolution, are taken into 
account. The GRB rate from our best result gives a total number of 4568^29 
GRBs per year that are beamed toward us in the whole universe. 


Introduction 


Gamma-ray bursts (GRBs) are one of the most energetic phenomena in the universe. 
Observationally, GRBs are usually categorized as long and short bursts, with an empirical 
separation of two seconds in their observational durations. The burst duration is often 
referred as T 90 , which is the light curve period that encloses 90% of the burst photons. Long 
GRBs are expected to result from explosions of massive stars with powerful central engines 


such as a black holes (e.g., Heger et al. 2003). Additionally, observations have shown that 


at least some long GRBs are connected to Type Ic supernovae (e.g., 

Galama et al. 

1998 

Bloom et al. 2002 Della Valle et al. 2003; Stanek et al. 2003 Thomsen et al. 2004; Campana 

et al. 2006 Chornock et al. 2010 Berger et al. 2011 ; [Vergani et al. 2011 Sparre et al. 2011 

Melandri et al. 2012 

). Short GRBs are believed to originate from compact-object mergers due 


to their different host galaxy properties and non-detections of the accompanied supernovae 
(e.g., Eichler et al.||1989 Nakar||2007 Fong et ahj|2010 , 2013). Here we will focus on the long 
GRBs, since they are closely related to massive stars, and hence trace the star-formation 
history more directly. Throughout the paper, GRBs refer to long bursts unless otherwise 
specified. 

Due to their extraordinary luminosities, GRBs provide a unique and independent method 
for measuring the cosmic star-formation rate (SFR), especially at large redshifts where it 
becomes difficult for other methods. Many efforts have been done to map out the intrinsic 
cosmic GRB rate as a function of redshift from current observations (e.g., Guetta & Della 


Valle 2007 Guetta & Piran 2007 Yiiksel et al. 2008 Kistler et al 

2008b Butler et al. 

2010 

Robertson & Ellis 2012; Pclangeon et al. 2008 Salvaterra et al. 

2009 Campisi et al. 

2010 

Wanderman & Piran|2010 Virgili et al. 1 201 1[ Qin et al. |2010 Salvaterra et al. 2012 Goward 


et al.||2013 Kanaan fe de Freitas Pacheco||2013). Results from these studies suggest that the 
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cosmic GRB rate generally follows the shape of the cosmic SFR at low redshift. However, 
at large redshift (z > 5), several groups have suggested a higher GRB rate than previously 


thought based on measurements of the SFR. (e.g., Le & Dermer 2007 Yiiksel et ah 2008 


Kistler et al. 2009 Butler et al. 

2010 

Ishida et al. 

2011 Tanvir et al. 

2012 

Jakobsson et al. 

2012 

). For example, 

Kistler et al. 

(2009 

) conclude that at high redshift the GRB rate does 


not trace the commonly adopted SFR from Hopkins & Beacorn (2006). These authors further 


state that the higher SFR implied by the GRB rate can be explained when including the star 
formation from undetectable galaxies at the faint end of the UV luminosity function, and 


therefore stars alone are sufficient to explain the reionization in the early universe (Kistler 


et al. 2013). Tanvir et al. (2012) used the high-redshift GRBs detected by Swift for locating 
their host galaxies and performing observations via the Hubble Ultra-Deep Field. Based on 
the non-detection of the host galaxies, they put an upper limit on the SFR of these galaxies 


and came to a similar conclusion as Kistler et al. (2009), that most of the star formation at 


high redshift comes from low-luminosity galaxies. 

To estimate the intrinsic GRB rate from current observations, one needs to convert 
the observed rate back to the intrinsic rate based on GRB luminosity distribution, the 
survey sensitivity, and other GRB characteristics such as the GRB spectral information, 
burst shapes/durations, and the beaming angles. Recently, several groups have done some 
careful examinations of the cosmic GRB rate via Monte Carlo approaches with adjustable 


parameters in the GRB rate and luminosity distribution for fitting (Butler et al. 2010 


Wander man & Piran 

2010; 

Qin et al. 

2010 

Virgili et al. 2011 


Some studies also adopt 
additional correlation between the GRB characteristics (e.g., luminosities, the peak energy 


-^peak °f burst spectrum) in their Monte Carlo search (e.g., Butler et al. 2010). 


Most of these studies focus on searching for the intrinsic GRB rate that are beamed 
toward us, because we can only observe GRBs that are pointed at us, and it is extremely 
difficult to measure the beaming angle. Despite all the effort, there remain large uncertainties 
in the beaming angle and the beaming factor, which is the fraction of GRBs that are beamed 
toward us. The beaming factor can range from ~ 50 to ~ 500 (e.g., Frail et al. 2001 Guetta 


et al. 2005). To avoid involving further uncertainties in our calculation, we also focus on 
finding the intrinsic GRB rate for GRBs that are beamed at us. All the GRB rates in this 
paper refer to GRBs that pointed toward us, unless specifically noted. 


Most recent studies concerning GRB rate must estimate a sensitivity for the survey 
of each detector under consideration. A single detection threshold is most commonly used 
for estimating the survey sensitivity, in which a flux (or photon count rate) limit is used 
as the instrument limit and a GRB with flux (or photon count rate) above the limit is 
considered as a trigger event. This is generally a good assumption for GRB instruments 
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prior to Swift. Unlike previously flown GRB instruments, however, Swift adopts a much 


more complex trigger algorithm in order to maximize the number of GRB detections (Band 


2006). The Burst Alert Telescope (BAT) of Swift adopts over 500 “rate trigger” criteria, 
based on photon count rates, and additional “ image threshold” based on the real image 


generated for further confirmation and localization (Barthelmy et ah 2005 Fenimore et al. 


2003, 2004 McLean et ah 2004 Palmer et ah 2004). Therefore, a single-detection threshold 


approximation might not be proper for correctly estimating the survey sensitivity of Swift. 
A few studies thus adopt more sophisticated approximations of the complex trigger 


algorithm of Swift. For example, Butler et al. (2010) used an empirical combination of photon 


counts and pulse durations to approximate the signal-to-noise ratio cut of Swift detections. 


Wanderman & Piran (2010) introduced an empirical probability function of detectability for 
bursts with flux close to the assumed single flux threshold. Additionally, these authors also 
adopted a empirical function to account for the probability of redshift measurements. 


Despite all the deliberate methods adopted to explore the intrinsic GRB rate, it remains 
difficult to quantify the selection biases from the complex BAT trigger algorithm. We hence 
proceed with an alternative approach by actually simulating the BAT trigger algorithm, 
including all the rate trigger criteria and image threshold. We will use this “BAT-trigger 
simulator” to search for cosmic GRB rate and luminosity distributions that generate a mock- 
triggered sample that can reproduce the observational GRB characteristics. 


This paper is organized as follow: Section [2] describes the complex trigger algorithm 
adopted by the BAT. Section [3] explains the method we use in our simulations, including 
generating GRB light curves based on input burst properties and simulating the BAT trigger 
algorithm. Section [4] summarizes the observational GRB samples we adopted to be compared 
with our simulations. Section [5] presents our best finding of the cosmic GRB rate and 
luminosity functions that generate results which match the best with the observational GRB 
characteristics. Section [6] explores the consequences of adopting different distributions of 
GRB spectra (-Upg C ak ), including spectral evolution. Section [T] discusses the advantages of 
using the complex BAT trigger algorithm, and the selection biases introduced by using a 
single trigger criterion. Section [8] shows the GRB detection fraction as a function of redshift, 
estimated by the BAT-trigger simulator. Section [9] compares our best-fit cosmic GRB rate 
with the cosmic SFR. Section [TO] explores the possibility of luminosity evolution. Finally, the 
results and their implications are presented in Sect. [lTj Throughout this study, we adopt 


a standard flat ACDM model with f2 m = 0.274, = 0.726, and H 0 = 70.5 km s 1 Mpc 

(Spergcl et al. 2007). 


,-i 
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2. Swift- BAT’s trigger algorithm 


During the observational process, BAT constantly calculates the signal-to-noise ratio 
of the detected light curves based on photon count rates in some assigned foreground and 
background time periods. Each foreground and background period are separated by an 
“elapse time” period to reduce the chance of the background photons being contaminated 
by the foreground photons. An event is “triggere” when the signal-to-noise ratio calculated 
within a given interval exceeds the assigned signal-to-noise ratio threshold. We adopt Eq. 3 in 


Graziani (2003) for calculating the signal-to-noise ratio. This is generally the same equation 


used by the BAT flight software, except the equation in BAT’s algorithm includes an extra 
term to prevent errors occurring when there are zeros in the denominator, that is, zero 
photons in the background light curves (Fcnimore et ah 2003). In our simulation we will not 
have the case with zero photons in the light curves, hence we do not need to include such 
term. 


BAT has 674 different rate trigger criteria. Each rate trigger criterion adopts a different 
signal-to-noise ratio threshold, different foreground, background, and elapse time periods for 
calculating the signal-to-noise ratio, and covers a different energy band (Fenimore et al. 2003). 
The complex trigger algorithm is implemented to increase the chance of correctly bracketing 
the GRB pulse shape and hence successfully find a GRB. After an event is triggered by one of 
the rate-trigger criteria, an image will be generated for further confirmation and localization. 
During this imaging process, an additional signal-to-noise ratio based on the real image will 
be calculated. The rate-triggered event will be confirmed as a real detection if (1) the image 
signal-to-noise ratio is higher than the image threshold (signal-to-noise ratio > 6.5), and 
(2) the event is compared with current on-board sky catalog and no known source is found 
to be at the event location. Out of all the rate trigger criteria, there are 180 criteria that 
have extremely short foreground periods (< 0.064 sec). These short trigger criteria are 
particularly aimed for detecting GRBs with very short durations. Since we focus on long 
GRBs in this paper, we do not include these short trigger criteria in our simulations. 


A throttling process is implemented to reduce the number of trigger criteria if the CPU 
is getting too far behind in processing the data. Additionally, BAT will also temporarily 
suspend the trigger process during some particular occasions, such as satellite slewing, when 
there are too many hot pixels that significantly reduce the detection sensitivity, or when 
the satellite is passing through the South Atlantic Anomaly region, which is an area of 
the Earth’s orbit that contains much higher background flux. In our simulation, we only 
consider the trigger process during normal situations. However, we take into account these 
“deadtimes” of observations by multiplying the simulated trigger rate by a fraction of the 
survey time (~ 90%) when the telescope is actually performing the observations. 
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In addition to the rate trigger process, a burst could also be found by an independent 
“image trigger” process, in which a regular image is generated by BAT every minute or 
longer to search for dim GRBs that are missed by the rate trigger. A typical image trigger 
process generates images from light curve periods of one and five minutes. However, when 
a rate trigger is active, the flight software also generates images from extended periods with 
increments of 8 seconds. That is, images from light curve periods of (64 + multiples of 8) 
seconds (i.e. , 72s, 80s, 88s... etc) will also be created, until the extension time reaches another 
minute. 


3. Simulations of GRB light curves and BAT trigger algorithm 

In order to search for a more robust cosmic GRB rate and to study possible systematic 
effects due to the complex trigger algorithm of BAT, we developed a code that is capable 
of creating mock observed GRB light curves based on adjustable burst properties, and sim- 
ulating the BAT trigger algorithm for the first time, including simulating hundreds of rate 
trigger criteria and mimicking the image threshold. This program contains three main parts, 
as described in the following subsections. 


3.1. Creating the mock light curves in the GRB rest frame based on assumed 

GRB characteristics 


We create a sample of mock GRB light curves in the rest frame based on several input 
GRB properties, which include: 


1. Cosmic GRB rate. We adopt the functional form in Wanderman & Piran (2010), 


Rgrb(z) — Rgrb{z — 0 ) 


( 1 ) 


which assumes a simple broken power law shape that generally follows the shape of 
the cosmic SFR, with adjustable parameters. That is, the cosmic GRB rate increases 
to some redshift zi, and decreases afterward, as shown in Eq. [TJ 

(l + z) n \ z<z u 

(1 + z 1 ) ni ~ n2 (1 + z) n2 , Z>Z 1 

Rgrb(z) is the comoving rate dN/ dV comov dt SIC , where dV coraov is the comoving volume 
and dt src is the time interval in the source frame. -Rgrb(^) can be converted to the 
observed rate i?GRB;dz = dN / dfldzdt Q b s in unit of number per solid angle (dfi), per 
redshift interval (dz), and per time interval in the observed frame {dt Q b s ) by 

dVr.omovdt src RgRb(z) dV c omov 


-RGRB;dz(^) ~ -RgRbW 


dfldzdtobs 


(1 + z) dQdz 


( 2 ) 


7 


2. GRB luminosity function: Again, we adopt the functional form in Wanderman & Piran 


(2010), which is also a simple broken power law function, 


m 


dN 

dL 


(£)*, l<u, 
(£) w , L>L* 


( 3 ) 


where x, y, and L* are adjustable parameters. Note that the luminosity L refers to the 
peak luminosity, not the average luminosity. 


3. GRB pulse shape: Using the observed GRB light curves from the real BAT-detected 
GRBs with known redshifts, we creat a library of rest-frame GRB light curves that 
contains 139 different GRB pulse shapes. Data with signal-to-noise ratio higher than 
3cr are considered as part of the GRB pulses, while data with signal-to-noise ratio 
below 3a are considered as noise and thus are ignored when constructing the GRB 
pulse shapes. The light curves range from duration T 90 = 2.16 sec to 658.2 sec. We 
randomly choose the GRB pulse shape from this library to create simulated light curves 
in the rest frame. 


3.2. Simulating the observed light curves with accurate energy response of 

the BAT 


We convert the GRB light curve in the rest frame created in part 1 (Sect. 3.1) into the 
observational light curve in units of photon count rate. The conversion is calculated using 
XSPECQ a program that is capable of converting flux into photon count rate based on the 


energy response of the instrument with an input GRB spectrum ( 

Arnaud 

1996 

). We adopt 

the commonly used Band function as the GRB spectrum (Band et al. 

1993 

), which has the 


form 


dN 

~dE 


KiE« exp[-E(2 + a)/EgU, 

K-2 E\ 


E < 


(“- 0 ) 


peak 


E > 


2+a ’ 


( 4 ) 


— 2+a 


Epea k i s the peak energy in the vF v spectrum in the observer frame. The normalization 
factors K\ and K 2 is decided by the GRB luminosity after redshifting into the observer 
frame. We assign different a and (3 in our mock GRB sample based on the observed spectral 
distribution reported in 


Sakamoto et al. 


(2009). We leave the E°^l k distribution to be 


one 


of the adjustable functions to explore possible consequences when different distributions are 
used. 


1 http: / /heasarc. nasa.gov / xanadu / xspec / 


We also take into account different energy responses of bursts coining from different 
angles relative to the detector plane. In other words, for the same burst, the simulation will 
create a light curve with higher photon count rate if the burst is on-axis, and lower photon 
count rate if the burst is off-axis. Moreover, BAT separates the detector into four quadrants 
and can trigger an event based on the light curve recorded in only some quadrants of the 
detector. The reason for doing this is to reduce noise for those GRBs coming from off-axis 
angles that illuminate only part of the detector. By calculating the signal-to-noise ratio of 
only the illuminated part of the detector, BAT can get a higher signal-to-noise ratio for these 
bursts and hence increase the chance of triggering these off-axis bursts. 

Table 1: Relation of the grid ID and the burst incoming angle (in degrees) relative to the 
detector’s plane, with zero degrees indicating an on-axis burst. There are no grid ID 6 and 
28. 


Grid ID 

Angle (°) 

Grid ID 

Angle (°) 

Grid ID 

Angle (°) 

Grid ID 

Angle (°) 

1 

50.69 

10 

19.27 

18 

26.56 

26 

46.63 

2 

40.68 

11 

31.39 

19 

44.99 

27 

56.99 

3 

35.02 

12 

46.63 

20 

56.29 

29 

50.66 

4 

40.68 

13 

56.99 

21 

56.99 

30 

40.68 

5 

50.66 

14 

56.29 

22 

46.63 

31 

35.02 

7 

56.99 

15 

44.99 

23 

31.39 

32 

40.68 

8 

46.63 

16 

26.56 

24 

19.27 

33 

50.66 

9 

31.39 

17 

0.076 

25 

31.39 




We label the BAT detector plan with some grid ID numbers. GRBs with different 
incident angles will fall on different locations on the detector plane with different “grid ID”, 
and different calibrations are used based on the “grid ID” of the burst. Figure [I] plots the 
relative location of each grid ID on the detector plane. Table [l] shows how the incident 
angle corresponds to each grid ID. An incoming angle of zero degrees corresponds to an 
on-axis burst, that is, the burst is directly above the detector plane. A larger incoming angle 
indicates that the burst is more off-axis. Grid ID 6 and 28 do not exist due to historical 
naming reasons. To simulate the actual trigger algorithm, our program also simulates the 
partial illumination factor based on the “grid ID” (i.e. , incident angle) of the burst, and 
creates light curves for four different quadrants of the detector accordingly. Figure [2] shows 
an example of simulated GRB light curves of the same burst with different grid IDs. One can 
see from the figure that the incident angle has a significant effect on the detector sensitivity. 
An off-axis burst can be ~ 20 times dimmer than an on-axis burst, and hence become 
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Fig. 1. — Location of the Grid ID on the detector plane. X and Y are the detector coordinate 
with arbitrary unit. The center of the detector is located at X=0, Y=0. 
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Fig. 2. — Simulated GRB light curves of the same burst with different incident an- 
gle relative to the detector plane. Specifically, bursts in the figure with incident angles 
= {0°, 30°, 45°, 55°} are generate from bursts with grid ID = {17,16,15,14}, respectively. 
Light curves are binned in 0.64 second. 
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undetectable. 


In addition, the active number of BAT’s detector changes with time. BAT has a total 
number of 32768 detectors. However, during the eight years of operation, many detectors 
are getting noisier and are automatically or manually turned off. Therefore, the average 
number of active detectors decreases each year, as shown in Fig. [3} This factor is taken into 


32000 r 


130000 


a! 28000[ 

Q 260001 
01 

~ 240001 
u 

^ 220001 


53 20000[ 

_Q 

^ 180001 


16000 


2005 2006 2007 


2008 2009 

Year 


2010 2011 2012 


Fig. 3. — The average number of BAT’s active detector as a function of year. 



Fig. 4. — Simulated GRB light curves of the same burst with different number of active 
detectors. Light curves are binned in 0.64 second. 
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account in our program by simulating light curves with different energy responses according 
to different number of active detectors. A burst would have a simulated light curve with 
higher photon count rates if we set a larger number of active detectors. Fig. [4] shows an 
example of how the simulated light curve changes with different number of active detectors. 


Furthermore, many observed GRBs show time-evolving spectra. Therefore, we imple- 
ment the option to include spectral evolution in our simulation. A majority of well-studied 
GRBs show a hard-to-soft spectral evolution as bursts become dimmer (e.g., Liang & Kar- 


gatis |1996 Crider et al. 1999 |Ryde fe Svensson 2000 |Zhang et al. |2007 Racusin et al. 


Yonetoku et al. 

2008; 

Page et al. 

2009; 

Filgas et al. 

2011) 


Some 


studies proposed that the spectral evolution is related to the pulse shape (Liang & Kargatis 


1996 Crider et al.||1999 Ryde fe Svensson||2000 ) . However, there are also bursts that do not 


show strong spectral evolution nor a monotonic behavior (Zhang et al. 2007). For simplicity, 
we adopt the assumption that evolves with the pulse shape. In other words, 

increases as a GRB becomes brighter, and decreases as the burst fades. Fig. [5] demonstrates 
an example of adding A°g® k evolution in the simulated light curve. 


Finally, we add some background noise to the simulated GRB light curves. We create 
a library of 371 background levels from the light curves of real BAT-detected GRBs. The 
backgrounds added to the simulations are randomly chosen from this library and are fluc- 
tuated with Gaussian noise. The lower four panels of Fig. [6] show several examples of the 
simulated light curves with background noise in the observational frame with units of photon 
count rate. The light curves shown here are in the 25-100 keV energy band, and are the 
summaztions of all four detector quadrants. Each panel shows the same burst with differ- 
ent incident angles, corresponding to the original light curve in Fig. [2] (without background 
noise). The example shows that the 55° off-axis burst is now completely buried under the 
background noise, and thus is undetectable. The grey and blue shaded regions indicate the 
foreground and background periods from the trigger criterion that detects the burst in our 
simulation (see the following section, Sect |3.3.1| for more discussion about the trigger simu- 
lator). For comparison, the top panel of Fig. [6] shows an example of a real GRB light curve 
(GRB120701A). The real GRB light curve is also a summation of light curves from all four 
detector quadrants, and in the energy band 25-100 keV. Note that real GRB light curves 
might not have flat backgrounds like those in our simulations, as shown in this example. 
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Fig. 5. — An example of the simulated GRB light curve with spectral (E^^) evolution. 
Light curves are binned in 0.64 second. 

3.3. BAT-trigger simulator that simulates the complex trigger algorithm of 

the BAT 

The code we developed is capable of simulating both the rate trigger and image trigger 
processes. The following sections describe in detail the methods we adopt to simulate the 
BAT-trigger algorithm. 
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Fig. 6. — An example of the simulated GRB light curve with background noise. For 
comparison, a real GRB light curve is also shown in the top panel. Simulated light curves 
with different incident angles are plotted in different panels. These are the same light curves 
shown in Fig. [2j but now with background noise added. The grey and blue shaded areas on 
top of the simulated light curves indicate the foreground and background periods that give 


the highest signal-to-noise ratio using our trigger simulator (see Sect. 3.3.1). Light curves 
are binned in 0.64 second. 
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3.3.1. Simulating the first part of the trigger process: the rate trigger 


Our program follows the same procedure adopted by BAT for the rate trigger. That is, 
we move through the light curve and calculate the signal-to- noise ratio of each bracketed time 
period (rb ra cket) specified by a rate trigger criterion. Each bracketed time period contains one 
foreground period, one or two background periods (depending on the trigger criterion), and 
an elapse time period between the foreground and background periods. For those trigger 
criteria that contain two background periods, the background periods Tbi and are placed 
before and after the foreground period r f , respectively. Both of the background periods are 
separated from the foreground period by the elapse period r e , to make sure the background 
light curve does not include photon counts from the source. For the trigger criteria that only 
use one background period Tbi, it is placed before the foreground period. An elapse period 
r e is also placed between Tbi and Tf for these one-background criteria. 


We follow Graziani (2003) in defining the signal-to-noise ratio, 


/if 


Tbi + 7b2 2 

Tf 



Signal-to-Noise Ratio 


72f 

(hf + S 2 )V 2 ' 


(5) 

( 6 ) 
(7) 


The signal-to-noise ratio is calculated by Eq. [?} using all the time periods described above 
and the corresponding photon counts in those periods. In the equations, nj,i, and rif 
represent the photon counts in Tbi, Tb 2 , and Tf, respectively. For the criteria with only one 
background, the signal-to-noise ratio is calculated by this same equation, with all the terms 
related to the second background (i.e. , rb 2 , dropped. 


After the signal-to-noise ratio is calculated, the program moves to the next time step, 
which in our simulation is set to be the next light curve bin, and calculates the signal-to-noise 
ratio in that time period. The program moves through the whole light curve and records 
the time periods where the signal-to-noise ratios are higher than the threshold determined 
by the specific trigger criterion. 


Each rate trigger criterion has different Tbi, r b 2 , Tf, r e , and signal-to-noise ratio threshold 
for triggering. Additionally, each trigger criterion uses light curves from different energy 
bands, and light curves from different summations of the four detector quadrants. The 
program thus goes through the light curves multiple times with settings specified by each 
criterion. The trigger criteria we adopt are identical to those used by the BAT flight software. 
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3.3.2. Mimicking the second part of trigger process: the image threshold 

After an event is triggered by one or more rate trigger criteria, BAT creates real sky 
images for further confirmation and localization, ffowever, creating images is very time- 
consuming and requires substantial computational power. Fortunately, it is possible to 
mimic the image threshold without creating real images in our simulations. Instead, we use 
information from the light curves and results of the rate triggers. This method allows us to 
calculate the image threshold without high demands of computational resources, and thus 
to be able to explore a larger parameter space in our search of the GRB characteristics. 

This approximation is feasible in our simulation because of the following reasons. For 
a real burst observation, images are created to obtain better localization information and 
clarify that the event is not from a known stellar object. However, in our simulation all the 
sources generated are GRBs. Hence, there is no confusion with other sources and all we 
need to know is whether the burst is triggered or not. Another reason for BAT to create a 
real image is to determine the real background level at the time when the burst happens. In 
reality, the background is constantly changing with time. Therefore, the background level 
when the burst happens is likely to be different than the background level before or after the 
burst. However, we assume a flat background with Gaussian noise in our simulation, and 
thus the background level should be time independent and the photon count rate calculated 
from the background period in the rate trigger should be the same as that obtain by an real 
image. Therefore, we can mimic the image threshold based on information gathered from 
the rate trigger. 

When a rate trigger criterion successfully triggers an event, BAT uses the accumulated 
photon counts in the foreground period of that rate trigger criterion to create an image. The 
photon counts used to create this image are added from all four quadrants of the detector, 
even if the rate trigger criterion that triggered this event only uses light curves from part of 
the detector. Therefore, when mimicking the image threshold, we use the foreground and 
background periods from the successful rate trigger criterion, but use the light curves from 
all four quadrants to calculate the signal-to-noise ratio. 

Determining the corresponding signal-to- noise ratio of the image threshold using only 
information from the light curve can be a little bit complicated. In reality, the image signal- 
to-noise ratio is calculated from the real image by fitting a point spread function. Therefore, 
the signal-to-noise ratio might not be identical to the signal-to-noise ratio calculated using 
light curves by Eq. [7j In an ideal case of flat background and only one source (i.e., the burst) 
per image, the signal-to-noise ratio calculated from point spread function fitting should have 
a one-to-one correlation to the signal-to-noise ratio calculated from the light curve. Thus, we 
can approximate the image threshold if we can find the corresponding signal-to-noise ratio 
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calculated from a light curve to that estimated from an image. That is, we need to know 
how to convert the image signal-to-noise threshold to a threshold in the light curve domain. 

To determine the corresponding signal-to-noise ratio of the image threshold, we apply the 
BAT-trigger simulator to 121 light curves of real BAT-detected GRBs. Using the foreground 
and background periods of the successful trigger criteria, we calculate the signal-to-noise ratio 
using light curves cladded from all four quadrants. Independently, we create real images for 
these real GRBs using the foreground period determined by the successful trigger criteria 
and calculate the image signal-to-noise ratio using the ground software of BAT (IIEASoftQ). 
We plot the signal-to-noise ratios calculated by these two methods to End the correlation in 
Fig.0 As discussed earlier, the scatter of the points shown in Fig. [7] is likely due to the large 



Fig. 7. — Correlation of SNR image , the image signal-to-noise ratio using the BAT software 
(HEASoft), and SNR 4quad , the signal-to-noise ratio from the four quadrant light curves. Red 
line shows the y 2 fitting, which has the form of SNR image = 2.47 + 0.44 SNR 4qua d. 

fluctuation of the backgrounds of these real GRB light curves, and/or from contamination of 
other sources in the light curves and images. For a first-order approximation, we simply fit a 
linear function to the correlation of the image signal-to-noise ratio calculated from the BAT 


2 http://heas arc. nasa.gov/lheasoft/ 
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software HEASoft (SNR image ) and the signal-to-noise ratio calculated from photon count rate 
using all four quadrant light curves (SNR 4quad ). This fitted line has the functional form of 
SNR image = 2.47 + 0.44 SNR 4qua d, with reduced y 2 = 11-12 (degree-of-freedom = 118). The 
signal-to-noise ratio for the real image threshold is ~ 6.5 to 7.0, which corresponds to a 
signal-to-noise ratio of ~ 10 calculated from light curves. Thus, we adopt the signal-to-noise 
ratio = 10 to be our threshold for mimicking the image threshold. 


3.3.3. The image trigger 

In addition to the regular two-stage trigger algorithm (rate trigger followed by image 
threshold), BAT also creates images regularly in the 15 — 50 keV energy range for longer 
durations (> 1 minute) to search for bursts missed by the rate triggers. GRBs found in this 
independent image trigger process usually have low fluxes but high fluences, which are likely 
due to relatively long and slow-changing light curves. These GRBs are thus hard to detect 
using rate triggers with shorter trigger durations. 

The signal-to-noise threshold for this image trigger process is ~ 7.0 to 7.5, which is 
similar to the image threshold following the rate trigger. Therefore, we adopt the same 
method as discussed in the previous section with the same criterion of signal-to-noise ratio 
~ 10 for mimicking the image trigger. However, we modify the foreground and background 
durations to replicate the longer exposure time in each image. 

As mentioned in Sect. [2j the BAT flight software makes images with many different 
durations, with some durations only available when a rate trigger is active, and thus gives 
some randomness to the ranges of exposure time of these real images. In order to determine 
the foreground durations for the image trigger process in our simulation, we list the time 
intervals of the image exposure times for the real BAT image-triggered GRBs, and adopt 
them as the foreground durations in our simulation. These durations are 64, 72, 88, 120, 
128, 168, 192, and 320 seconds. The background durations are set to be the same as the 
foreground durations in our simulation to calculate the signal-to-noise ratio of the bursts. 
We put in a 32 second elapse time between the background and foreground periods, to make 
sure the background is not contaminated by the burst light curve when there is a detection 
in the simulation. We only consider triggers in the 15 — 50 keV band for the image trigger, 
because this is the only energy band used by BAT during this process. 
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3.3.4 ■ Comparing our simulation with BAT’s sensitivity 


To test whether our program correctly simulates the complex BAT-trigger algorithm, we 
compare the GRB peak fluxes of the “triggered” bursts in our simulation to those measured 
from the real GRBs detected by BAT. 


Panel (a) of Fig. [8] shows the peak fluxes of real BAT-detected GRBs with respect to 
the grid ID of the detector plane. As described in Sect |3.2 the grid IDs are simply the 
number labels on the detector plane, and thus correspond to the incoming angles of the 
bursts relative to the normal axis of the detector. Each point in the figure represents one 
burst. The fluxes of these bursts are adopted from Sakamoto et al. (2011). Blue dots in the 
plot indicate GRBs detected by rate triggers, while red crosses represent bursts found by 
image triggers. 


Panels (b) and (c) of Fig. [8] plot the fluxes from the mock “triggered” bursts in our 
simulations with 20000 and 30000 active detectors, respectively. Again, blue dots show the 
rate-triggered bursts, while red crosses indicate the image-triggered events. Because the 
number of active detectors decreases with time, the sensitivities using two extreme numbers 
of active detectors are plotted for comparison. Results show that decreasing the number 
of active detectors from 30000 to 20000 reduces the sensitivity by a factor of ~ 3, and has 
more significant impact on off-axis bursts than on-axis events. Moreover, the sensitivity of 
the image trigger is less affected by reducing the number of active detectors. The input 
GRB characteristics of this mock sample are based on the best result in our search, which 
is summarized in Table [2j Sect. [5] This sample creates plenty of bursts that have fluxes 
in the range of ICR 9 erg s” 1 cm -2 to ICR 6 erg s _1 cm -2 . Therefore, the non-detection of 
low-flux bursts is due to the sensitivity of our trigger-simulation program, instead of the lack 
of low-flux events. 


Results show that the rate trigger process in our trigger simulator can detect GRBs 
with fluxes ~ ICR 8 erg s -1 cm -2 for directly on-axis bursts, and fluxes ~ ICR' erg s^ 1 enr 2 
for extremely off-axis events, which is very similar to the real BAT sensitivity, ft is harder 
to compare the sensitivity of the image trigger process, due to the very low number of 
statistics in real BAT detections. However, in general the image trigger process in our 
simulations detects bursts with fluxes that are slightly lower than those found by the rate 
trigger algorithm, as expected. A similar trend is also seen in the real BAT-detected GRBs. 
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Fig. 8. — Panel (a): 1-s peak energy flux vs. grid ID for the BAT detected bursts from 
GRB041223 to GRB091221 (Sakamoto et al. 2011). Panel (b): Peak energy flux vs. grid 
ID for the triggered GRBs in the mock sample, assuming 20000 active detectors. Panel (c): 
Peak energy flux vs. grid ID for the triggered GRBs in the mock sample, assuming 30000 
active detectors. 


4. Observational distributions of the GRB characteristics 


In our search for the intrinsic GRB characteristics, we need to compare the GRB proper- 
ties from our mock-triggered sample, which contains the simulated bursts that are triggered 
by our trigger simulator, to the GRB properties from the real BAT-detected bursts. In 
this section, we describe the observational GRB samples we adopt to be compared with our 
simulated bursts. We also discuss the challenges and uncertainties in these observational 
measurements. 


The main GRB characteristics we use for comparison are the redshift and the peak-flux 
distributions. We iteratively adjust the input parameters until the simulated results match 
well with both the observed redshift and peak-flux distributions. The peak-flux distribution 
is chosen as one of the main guides in our search because the peak flux can be measured 
directly from observations and thus is less uncertain, as described in Sect. |4.2| The redshift 


distribution is also selected because our main goal here is to find the intrinsic GRB rate. 
In addition, we also consider E°^l k and E iso from real observations. However, due to the 
difficulties in the measurements of 77°g® k an d Also an d their large uncertainties in the BAT 
data, these two properties are only used as references but not as the primary guides in our 
search. 


In all these observed GRB samples we adopted, the bursts found by ground analysis, 
instead of flight software, are removed. This is because our trigger simulation follows the 
method of the flight software. Moreover, currently the ground analysis is mainly done by 
human and hence is less systematic. There are only a few ground-detected bursts in the 
observed GRB samples (zero in the redshift sample; 7 out of 409 in the flux sample; 9 out 
of 423 in the A pea k sample; see Sect 4.1 to Sect 4.4). Therefore, excluding these bursts does 
not have significant effect on the results. 


4.1. The redshift distribution of the BAT-detected GRBs 


Around 30% of all the BAT-detected GRBs have measured redshifts (Gehrels & Meszaros 


2012). These redshifts are usually measured from the burst afterglows and/or host galax- 
ies observed by follow-up ground observations. Therefore, the GRB redshift sample suffers 
greatly from observational biases, and is highly skewed toward lower redshifts. As there are 
many unpredictable factors that affect the redshift measurements, such as the local weather 
conditions and the number of telescopes that are available for follow-up observations, it is 
difficult to construct a complete redshift sample of the observed GRBs. 


Despite the complications, many studies have discussed the observational selection biases 
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and their effect on the GRB redshift distribution, and explored possibilities to construct 


a complete GRB sample (e.g., 

Coward et al. 

2008 

Fynbo et al. 

2009 

Jakobsson et al. 

2012 

Hjorth et al. 

2012 Salvaterra et al. 

2012 

Coward et al. 

2013 

). In particular, 

Fynbo 


et ah (2009) carefully consider possible observational constraints on performing follow-up 
observations and create a GRB sample that is less affected by the observational biases. 
These authors set up a number of criteria to select GRBs that have optimal conditions for 
follow-up observations. For example, bursts that are too close to the Sun or the Moon, or 
have high Galactic extinction, are removed from the sample, because it is harder to perform 
follow-up observations for these bursts and thus these bursts need to be extremely bright to 
have measured redshifts. 

There are 79 GRBs with redshift measurements (either from afterglows, host galaxies, or 


both) in the statistical GRB sample compiled by Fynbo et al. (2009) (Table 73 in their paper). 
These authors note that bursts with redshift measurements from the optical afterglows are 
not representative for all the Swift bursts in their statistical sample. The GRBs without 
optical spectroscopy in their sample are likely to suffer from dust obscuration other than 
the Galactic extinction, such as dust in the GRB host galaxies. Similar conclusion has 
also been obtained by The Optically Unbiased Gamma-Ray Burst Host (TOUGH) survey, 
which performs a systematic search for the GRB host galaxies and found that the host 
galaxies of GRBs with no optical and/or near infrared afterglows are significantly brighter 


and redder than those with optical and/or near infrared afterglows (Hjorth et ah 2012). 
Although it is possible to use host galaxy detections to recover these missing redshifts of 
those bursts without optical afterglows, this approach introduces other biases due to host 
galaxy brightness and the optical survey sensitivity. 


We therefore adopt the GRBs from the statistical sample constructed in Fynbo et al. 


(2009) that have redshift measurements either from only afterglows, or from both afterglows 
and host galaxies. We exclude bursts with redshift measurements only from host galaxies, 
because low redshift bursts have a higher probability of having detectable host galaxies. We 
also exclude GRBs with photometric redshifts due to their large uncertainties. Additionally, 
we remove one burst that does not have E iso information (see Sect. 4.4). There are 66 bursts 
in this redshift sample we adopt. 


4.2. The peak-flux distribution of the BAT-detected GRBs 

The GRB peak fluxes in the BAT energy range 15 — 150 keV can be measured directly 
from the BAT observations using the fewest assumptions about the burst characteristics. 
Therefore, the 15 — 150 keV peak flux is the least uncertain among all the GRB properties 


22 


we consider here. We adopt the 15 — 150 keV peak fluxes of the real BAT-detected GRBs 


reported in Sakamoto et al. (2011); a total number of 402 bursts are given. 


4.3. The A°j)® k distribution of the BAT-detected GRBs 


Not all of the BAT-detected GRBs have measured A°g® k . 


Therefore, we use the A obs 


estimator found in 


Sakamoto et al. 


peak 


(2009) to estimate the burst A° b ® k based on the power- 
law index (T) of the burst spectra when fitted by a simple power-law model. The power- law 
indices are reported in Sakamoto et al. (2011). 


The -E°gak estimator in 


Sakamoto et al. 


(2009) can only calculate A° b ® k in a limited range 
of power-law indices (1.3 < T < 2.3), which corresponds to A° b ® k values inside the detectable 
energy range of the BAT. In the cases where the power-law indices are out of the range, we 


use the power-law indices to estimate whether A° b ® k lies below or above the BAT detectable 
energy range. A small power-law index (T < 1.3) indicates that A° b ® k falls above the BAT 
detectable energy range (A° b ® k > 150 keV), while a large power-law index (T > 2.3) implies 
that -E°eak is lower than the BAT detectable energy range (A° b ® k < 15 keV). There are 414 
bursts in the Apeak sample we adopt, in which 26 bursts are below the BAT detectable energy 
range, and 80 bursts are above the BAT detectable energy range. 


4.4. The A<, 0 of the BAT-detected GRBs 


The total energy output (i.e. , the fluence) E iso of an observed GRB is difficult to measure. 
Both the burst spectrum and redshift are required to calculate E iso . However, most of the 


bursts are lacking redshift measurements (see Sect 4.1 ) and/or well-constrained A° b ® k , leading 


len A°ga k lies outside of the BAT energy 


to a poor characterization of the spectra, especially w’ 
range. Therefore, even with a burst that has an observed redshift, estimating Aj SO requires 
one to extrapolate the measured spectrum to an energy out of the detectable range, and 
hence introduces uncertainties. In addition, an observation might not capture the complete 
light curve of a burst due to the background noise. In order words, it is likely that BAT only 
detects the tip of the burst and thus underestimates the total energy output of the event. 


Butler et al. (2007, 2010) report a list of A; 


values of the BAT-detected bursts. They 
estimate the A iso in the T 90 burst duration and energy range 1 — 10 4 keV. Due to the difficulties 
in directly measuring A iso , these authors adopt a Bayesian approach to estimate the values, 
with a prior A° b ;( k distribution following results from the observations of GGAO/BATSE 
(Preece et al. 2000), the primary GRB instrument prior to Swift. Robertson & Ellis (2012) 
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compile a list of BAT-detected GRBs that have measured redshifts and E iso values. The 
majority of the E [so values of their list are from Butler et al. (2007, 2010) with a few more 


from Sakamoto et al. (2011). In addition, these authors calculate Ei SO 

(see 


for 29 new bursts 


Robertson & Ellis 


from their fluence and redshift values based on numerous references 
2012). The E iso estimations become more uncertain for bursts with E°^l k lying outside of the 


BAT detectable energy range, because it is hard to pinpoint the turnover of the spectrum. 
Therefore, for the E iso sample, we only consider bursts with E°^l k in the BAT energy range 
and adopt the corresponding Ei SO values from the list in Robertson & Ellis (2012). 


5. Searching for the intrinsic GRB characteristics 
5.1. General methods 


We modify the parameters for the intrinsic GRB redshift distribution and luminosity 
function in Eq. [l] and Eq. [3] to search for a set of the parameters that generates a mock- 
triggered burst sample that matches the best with the observed GRB characteristics. Since 
we follow the same functional forms for the redshift and luminosity equations as those in 
Wanderman & Piran (2010|, we start our search around numbers reported by those authors. 


In our search, we simulate 10000 bursts for each set of parameters in order to have enough 
mock-triggered bursts to reduce statistical fluctuations. Our code was run on computers with 
4 Intel quad-core Q9650 processors at 3 GHz on the Scientific Linux release 5.9 (Boron). 
Without including E° y ^ k evolution, our code takes around 5 hours to simulate 10000 bursts 
and to run them through the trigger simulator, for light curves with bin size = 1.6 seconds. 
If evolution is included, the simulation takes ~ 3 days for a sample of 10000 bursts 

with light curves binned into 1.6 seconds. If we decrease the light curve bin size (i.e. , increase 
the number of light curve bins for each burst), the time required for the trigger code to go 
through the whole light curve range increases as the total number of light curve bins. 


The equations for the redshift distribution and luminosity function we adopt contain 
seven parameters total (see Eq. [I] and Eq. [3]). Therefore, to run a complete Monte Carlo 
simulation and search through the full parameter space is highly demanding of computational 
power and time. For example, if we explore a range of ten values for each parameter, we will 
have 10' combinations for the parameter set, and hence the full Monte Carlo simulations 
with burst light curves binned to 1.6 seconds will take ~ 5700 years for the simulations to 
finish. One possibility might be adopting the Markov Chain Monte Carlo method. However, 
due to the complexity in the GRB observables we are taking into account (e.g., cosmic 
GRB rate, luminosity function, -E°g* k distribution), it is difficult to fold a good and efficient 
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algorithm that guarantees speeding up the process significantly and converging to the correct 
answer. Therefore, instead of searching the full parameter space, here we only try to find 
at least one possible set of parameters that generates a good match with the observed GRB 
characteristics. 


Among all the long rate trigger criteria of the BAT, 250 criteria have foreground periods 
longer than two seconds. Since we focus only on long bursts in this paper, most of the bursts 
that are triggered by criteria with foreground periods shorter than two seconds should also 
be triggered by criteria with foreground periods longer than two seconds. Therefore, to speed 
up the search process for the GRB properties, we only run through the 250 criteria that are 
longer than two seconds in our main search. There could be a few scenarios where the long 
GRBs are “long” by virtue of consisting of multiple “short” spikes that are well-separated 
(longer than the foreground period of the trigger criteria). For these bursts, it is possible 
that they could only be triggered by criteria with foreground durations shorter than two 
seconds, ffence, once we find a parameter set that matches well with observations, we will 
rerun the sample with all the long trigger criteria to check whether the results remain a good 
fit with observations. 


Additionally, our main searches are done without including E°^l k evolution, and with a 
fixed -Speak distribution, for the purpose of speeding up the search process as well. To decide 
what kind of Apeak distribution to use in our main search, we perform some test runs at 
the beginning using several very different functions for the E^ ak distribution (see Sect 


6.1 


for detailed discussion of the choices of functions). We find that in general, assuming some 
kind of intrinsic relation between E^ k and the burst energy output (e.g., luminosity, A iso ) 
seems to create a better match with observations. Therefore, in our main searches we adopt 


an Ap(f ak distribution that follows the Yonetoku relation of A®g C ak and L peak (Yonetoku et al. 
2004) (see more discussion in Sect. 6.1). Once we find a sample that matches well with the 


observations, we modify the E^ ak distribution using several very different functions to see 
whether this would cause significant changes in the results. Similarly, we also run the same 
set of parameters with A°g® k evolution to see how the results might change. If the effect of 
adopting a different A p g C ak distribution and/or A° aak evolution turns out to be significant, we 
adjust the parameters in the redshift and luminosity functions and repeat this process. 


Furthermore, as discussed in Sect. 3.3.4 the sensitivity of BAT, and hence the sensitivity 
of our trigger simulator, is slightly different for different numbers of active detectors. Ideally, 
to simulate the decrease in BAT’s sensitivity, we need to run our simulation for different 
numbers of active detectors and take an average of all the results, ffowever, this can easily 
increase the search time to a point that it becomes impractical. Therefore, we start our 
search with one number of active detectors. Usually we start with ~ 25000 active detectors, 
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which is the medium number of detectors for the eight years of BAT’s operation. Once 
we find a set of parameters that matches well with the real observations, we perform the 
simulations for several different numbers of active detectors and take the average. Minor 
adjustments of the parameters might be needed until this averaged result fits well with the 
observational GRB characteristics. 

To see how good the mock-triggered sample matches the observed GRB characteristics, 
we compare several burst properties of the mock-triggered sample with those of the real BAT- 
detected GRBs. As discussed in Sect. [|J we use the redshift and peak flux distributions as our 
main guides. For these two distributions, we perform the KS test to quantify how good the 
matches are between the distributions from the mock-triggered sample and those from the 
real BAT-detected GRBs. We modify the parameters in the redshift and luminosity functions 
until the KS test (with uncertainties) gives a significance level above 90%. Additionally, we 
use the -Ep^k distribution and the E^ k -E [so relation as references. We try to make these two 
distributions match as well as possible to the observed ones. However, due to the large and 
hard-to-quantify uncertainties in E°^l k and E iso (see discussion in Sect. |tJ) , we do not require 
them to match perfectly, and do not perform statistical tests on these two distributions. 


5.2. Results of the best-fit parameters 

The parameter set shown in Table [2] contains our best fit parameters for Eqjljand Eq. [3j 
This set of parameters generates mock-triggered bursts that have a redshift distribution 
and peak-flux distribution that match the best with those from the real observed GRBs. 
Table [2] also lists the KS test values for both the redshift distribution and the peak-flux 
distribution. Both the number “D” in the KS-test and the significance level of “D” (i.e. , the 
probability of “D” larger than the reported value) are given in the Table. The value “D” 
is the key parameter in the KS test, which indicates the maximum distance between the 
cumulative distribution of the two tested samples. Smaller “D” implies a better match of 
the two samples. All of the real BAT-detected bursts we adopted for comparison are GRBs 
before 2009. Therefore, we run the simulations five times (each time generates 10000 bursts) 
with a different average number of enabled detectors from year 2005 to 2009. The results 
of the simulations shown in the following figures are generated from the total 50000 bursts 
with different numbers of enabled detectors. 

Figure [9] shows the redshift distribution from this best-fit sample. Panels (a) and (b) of 
Fig. i give the input redshift and luminosity distributions, respectively, of all 50000 bursts 
created. Panel (c) of the figure shows the comparison of the redshift distribution between the 
mock-triggered bursts and real observations. The red bars in panel (c) show the normalized 
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Table 2: Summary of the set of parameters that generates results that match the best with 
the observed GRB characteristics. 


Rgrb(z = 0)[Gpc 3 yr 3 ] 

Zl 

n\ 

n 2 

0.84 

3.6 

2.07 

-0.70 


L*[erg s 1 

X 

y 

i^peak distribution 

2Q52.05 

-0.65 

-3.00 

Modified Yonetoku Relation (Eq. 

8) 


KS-test (D) for 
z distribution 

Significance for 
z distribution 

KS-test (D) for 
peak flux distribution 

Significance for 
peak flux distribution 

4.77 x 10~ 2 

99.79% 

3.091^04 x 10 -2 

85.59^;io% 


numbers of the simulated bursts that are triggered by our trigger-simulator. The blue dots 
in panel (c) show the normalized numbers of real BAT detections. The error bars along the 
y-axis are the statistical errors (i.e., square root of the number in each bin). The x-axis error 


bars simply represent the bin size. As mentioned in Sect. 4.1, we adopt the GRB redshift 


sample reported in Fynbo et al. (2009) for onr comparison 


Figure [lO] shows the peak-flux distribution from the best-fit sample. The top panel gives 
the input peak-flux distribution from the whole sample of 50000 bursts. The bottom panel 
shows the peak-flux distribution of the simulated bursts that are triggered by our trigger- 
simulator. Again, this mock-triggered sample is given as red bars and the distribution of the 
real BAT-detected GRBs is shown as blue dots. The peak fluxes of the real BAT-detected 
GRBs are reported in Sakamoto et al. (2011). These authors also listed the errors for each 
measured peak flux. Therefore, we run a quick Monte Carlo simulation to see how the 
match of the two distributions (i.e., the value from the KS test) changes if one allows the 
peak fluxes to change within their error bars. Results show that the uncertainty in the peak- 
flux measurements can change the “D” value in the KS test from 7.21 x 10~ 2 to 2.05 x 10“ 2 , 
which correspond to a significance level of 3.66% and 99.69%, respectively. These are the 
uncertainties we listed in Table [2} 


Figure 11 plots the E peak distribution for this best-fit sample and its comparison with 
the distribution from real BAT-detected GRBs. Panel (a) plots the input distribution 

of the 50000 simulated bursts. Panel (b) shows the comparison of the distribution 

between the mock-triggered bursts and the real BAT-detected GRBs. Similarly, this mock- 
triggered sample is plotted in red bars and the distribution of the real BAT-detected GRBs 
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is shown as blue dots. The _E°g® k values of the real GRBs are estimated from 


(2009 


2011 


fall out of t 

pobs 


Sakamoto et al. 


, as described in Sect. |4.3[ Because it is hard to estimate _E°g® k when the values 
re BAT energy range (~ 15 — 150 keV), we make two large bins for bursts with 
-^peak outside of the BAT range. Inside the BAT energy range, the data are binned into two 
smaller bins. 

Panel (c) of Fig. 


11 


shows the correlation. 


Red dots in this plot represent 
values from the mock-triggered sample. Dark red dots show the bursts with E ° values 
inside the BAT energy range, while light red dots are events with -E°g® k values outside of 
the BAT energy range. Blue crosses show the bursts from the real BAT detections. Due to 
the the large uncertainties of the -Epeak an d -Eiso values for real bursts, only events with the 
observed _E°g® k inside the BAT energy range are shown in the plot (see detailed discussion in 
Sect. [4]) . The E iso values in the figure are integrated over T 90 , both for the real GRBs and the 
simulated bursts. The intrinsic E iso values of the real BAT-detected bursts are unknown due 
to background noise. Therefore, we also restrict the E iso to the T 90 range for the simulated 
bursts to have a fair comparison with real observations. 

Results from our search suggest that a slightly modified Yonetoku relation 

-0.5 


£peak = 1-8x1 2.34 x 10" 5 x . 9 Lpeak . 

p eak i 10 52 erg g-1 


( 8 ) 


creates a better match for the -E pea k distribution, as shown in Panel (b) of Fig. 11 Everything 
inside the bracket of Eq. [8] is the original functional from Yonetoku et al. (2004]) , but we 
use peak luminosity E pea k instead of the isotropic luminosity L iso . In our modified Yonetoku 
relation, E^ ak is 1.8 times larger than that produced by the original Yonetoku relation. 
The main reason for us to consider this modified Yonetoku relation is because this relation 
generates higher E^ k bursts for the same E pea k, and thus increases the number of detections 
of higher E^y. bursts and matches the observations better. 

This best-fit sample predicts that Swift should detect ~ 96 bursts per year, which is in 
good agreement with the average number of ~ 95 GRBs per year from 2005 to 2009 based 


on real BAT observations (Sakamoto et al. 2011). The predicted detection rate for Swift 


(. Rswift ) is calculated based on the following equation, 

F^Swift 77-GRB;dz X /detect X FOV X t ; 


survey ? 


( 9 ) 


where 77.GRB ; dz is the observed GRB rate in units of number per redshift per solid angle 
per time in the observed frame (Eq [2]) , f detect is the detection rate, as shown in the third 
column of Table [3j FOV ~ 2 sr is the field-of-view of the BAT (Barthelmy et al. 2005), and 


^survey ~ 


90% is the fraction of the time that BAT spends on searching for GRBs. 


(a) Input Redshift Distribution 



(b) Input L peak Distribution 
(Before Trigger) 



ipeak [erg s 1 ] 


2 ^c) Redshift Distribution of the Mock-Triggered Bursts 
• • Swift GRBs (Fynbo et al. 2009) 



Redshift z 


Fig. 9. — Panel (a): the redshift distribution of all 50000 simulated bursts. Panel (b): the 
peak-luminosity distribution of all 50000 simulated bursts. Panel (c): the redshift distri- 
bution for the mock-triggered bursts (red bars), which are those bursts that are triggered 
by our trigger simulator. The redshift distribution of the real BAT-detected bursts is also 
plotted for comparison (blue dots; Fynbo et al. (2009)). Error bars along the y - axis show 
the statistical errors in each bin. Error bars along the :r-axis represent the bin sizes. 
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Input Peak Flux Distribution (Before Trigger) 



Peak Flux [erg/s/cm 2 ] 


Peak Flux Distribution of the Mock-Triaaered Bursts 



Peak Flux [erg/s/cm 2 ] 


Fig. 10. — Upper Panel the peak-flux distribution of all the 50000 simulated bursts. Bottom 
Panel, the peak-flux distribution for the mock-triggered bursts (red bars). The peak-flux 
distribution of the real BAT-detected bursts are also plotted for comparison (blue dots; 


Sakamoto et al. (2011)). Error bars along the y - axis show the statistical errors in each bin. 


Error bars along the x-axis represent the bin sizes. 
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(c) E iso vs. E;™ k 

of the Mock-Triggered Bursts 

■ ■ Mock triggered sample 

(15 keV < E ° ^ <150 keV) 
x x Swift GRBs (15 keV < E ^ < 150 keV; 
Butler et al. 2007, 2010; 

Sakamoto et al. 2011; 

Robertson et al. 2012) 

Mock triggered sample (all) 



Fig. 11. — Panel (a): the E^ ak distribution of all 50000 simulated bursts. Panel (b): the 
-^peak distribution for the mock-triggered bursts (red bars). The ££g® k distribution of the 
real BAT-detected bursts are plotted (blue bars; Sakamoto et al. (2011)). Error bars along 
the y-axis show the statistical errors in each bin. Error bars along the x-axis represent the 
bin sizes. Small and large bins are bursts wth E°g® k inside and outside of the BAT energy- 
range, respectively. Panel (c): E iso versus E^ ak . Dark red dots and blue crosses show bursts 
with in the BAT energy range for the mock-triggered sample and the real GRBs, 

respectively. The full mock-triggered sample is shown as light red dots. E- lso values in the 
plot are integrated over T 90 . 
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As discussed in Sect. 5.1, we rerun this best-fit sample with the complete set of long 
trigger criteria, including those with foreground periods shorter than two seconds, in order 
to make sure the results remain a good fit with observations. We use 26884 active detectors 
in this run, which is the average number of active detectors from 2005 to 2009. Results show 
that comparing with the same sample (with 26884 active detectors) using only the trigger 
criteria longer than two seconds, adopting the complete set of trigger criteria, change the 
detection rate from 13.73% to 14.01%. The KS test significance for the redshift distribution 
changes from 99.02% to 99.32%, and the KS test significance for the flux distribution changes 
from 81.80 +q 5;55% to 89.78lg7 4 g 0 %. All of these changes are significantly smaller than the 
statistical uncertainty. Therefore, results using our best-fit parameters remain good matches 
with observations when adopting the full set of long trigger criteria. 


6. T^peak distribution and evolution 


6.1. Aj% ak distribution 


The intrinsic E^ ak distribution remains uncertain and controversial. Many studies have 
suggested possible correlations between the total energy output of the burst and the peak 
energy of the vF u spectrum in the burst rest frame E^^. These relations often attempt 
to relate E^ ak and E iso , 
of the burst 


In order to see the consequences of adopting different E^ ak distributions, we test several 
functions that have significantly different shapes compare to those used in our best-fit sample 
(described in Sect.[5j). The distributions we test include: (1) A flat E^ k distribution in linear 
space. (2) A flat E^ ak distribution in logarithmic space. (3) A Gaussian E^ ak distribution 
in logarithmic space, with average = 300 keV and a — 1. (4) A special function 


luminosity 


et al. 2004 


the total energy fiuence of the burst, or E^ k and L iso , the total 


(Amati et al. 2002 Amati 2006; Ghirlanda et al. 2004 Yonetoku 


0(logio(E"' k )) oc 


logio(^U) if < 10 keV, 
log 10 (10 keV) otherwise. 


( 10 ) 


This function contains a lower number of bursts for 17pg C ak < 10 keV and follows a flat 
distribution in logarithmic space for events with E^ ak > 10 keV. 77°g ak evolution is not 
included in these simulations. 


Table [3] summarizes the major results from using different E^ ak distributions. All the 
results, except our best-fit sample (the one with a modified Yonetoku 77pg ak distribution), 
are based on 10000 simulated bursts and using 26884 active detectors, which is the average 
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number of active detectors from year 2005 to 2009 (the time period of the real BAT-detected 
bursts adopted in this paper). Our best- fit result contains a total number of 50000 simulated 


bursts with different numbers of enabled detectors, as described in Sect. 5.2 


l.Or 


(1) Flat i?° p b a y Distribution in Linear Space 

• • Swift GRBs (Sakamoto et al. 2011) 
Mock-Triggered GRBs 


l.Or 


(2) Flat Distribution in Log Space 

• • Swift GRBs (Sakamoto et al. 2011) 
Mock-Triggered GRBs 



(3) Gaussian E°^ v Distribution in Log Space 

• • Swift GRBs (Sakamoto et al. 2011) 
Mock-Triggered GRBs 


(4) Specific E °™ k Distribution (Eg. 10) 

• Swift GRBs (Sakamoto et al. 2011) 
Mock-Triggered GRBs 



SSc [keV] 


Fig. 12. — Comparison of the E°g® k distributions for the mock-triggered bursts (red bars) 
assuming different input E distributions (Sect. 6.1). 


The -E°g® k distribution of the real 


BAT-detected bursts are also plotted for comparison (blue bars; Sakamoto et al. (2011)). 
The narrow bins indicate bursts with E°^ k values inside the BAT detectable energy range. 
The two large bins contain bursts with E° g® k values outside of the BAT energy range. Error 
bars along the y-axis show the statistical errors in each bin. Error bars along the x-axis 
represent the bin sizes. 

For all the four different E^ ak distributions we test, the resulting redshift and peak flux 
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(3) Gaussian E^ k Distribution in Log Space 

• Mock triggered sample 
(15 keV < E°^ < 150 keV) 
x Swift GRBs (15 keV < E°^ k < 150 keV) 
Mock triggered sample (all) 
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(4) Specific E^!* v Distribution (Eg. 10) 

■ Mock triggered sample 
(15 keV < E°^ k <150 keV) 
x Swift GRBs (15 keV < E'^ k < 150 keV) 
Mock triggered sample (all) 
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Fig. 13. — Comparison of the resulting E iso versus E^ ak f° r the mock-triggered bursts using 
different input E^ aak distributions (Sect. 6.1 ). Dark red dots and blue crosses represent bursts 
with f?°gak hi the BAT energy range for the mock-triggered sample and the real BAT-detected 
GRBs, respectively. E°^l k for real GRBs are from 
(2011), and Robertson & Ellis (2012). For comparison, the 
also shown as light red dots. E iso in the plot represents E iso in the r 90 range. 


Butler et al. 

(2007, 

2010) 

Sakamoto et al. 


fill mock-triggered sample is 
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distributions of the mock-triggered samples remain good fits to the real observations. In 
fact, for the peak-flux distributions, for which we can quantify the uncertainties of the fit, 
the KS-tests show that all these matches based on different E^ k distributions are within 
the calculated uncertainties of each other. In other words, all the E^ k distributions in our 
tests give the same level of good fits to the observations, and thus comparing the results to 
the redshift and peak-flux distributions alone would not be sufficient to distinguish different 
-^peak distributions. However, using different E^f ak distributions does change the detection 
rates. For the four different cases we tried, the detection rate can change from ~ 4% to 
~ 14%, and thus will affect the normalization of the GRB rate (i.e. , the GRB rate at z — 0) 
up to 4 times higher. 


The major distinctions of adopting different E^ r e c ak distributions appear in the E peak 


12 


and 


13 


compare the resulting A°k ak 


distributions of the mock-triggered samples. Figure 
distributions and the E\ 

As one can see from these figures, all these additional Apeak distributions we test seem to 
generate worse matches to the observational sample than the modified Yonetoku sample 


Epe^k correlations based on these four different E^ ak distributions. 


see 


Fig. 11), in both the A° aak distributions and the A iso -Apg C ak correlations. However, due to 
the large and hard-to-quantify uncertainties in the A° aak and Ei SO of the real observational 
sample, it remains ambiguous whether any of these E^f ak distributions can be excluded. 
Therefore, these plots are only shown to indicate how the distributions can change if one 


assumes different Apg C ak distributions; no conclusion about the intrinsic Apg ak distribution 
can be drawn until direct measurements of E°^l k and A iso become available. 


Table 3: Summary of the results from simulations without A° aak evolution. The five samples 


shown here are based on different E^f ak , as described in Sect. 6.1 


zpsrc 

■^peak 

Detection 

Prediction for 

KS-test significance 

KS-test significance for 

Distribution 

Rate 

Swift [yr _1 ] 

for z distribution 

peak flux distribution 

Modified Yonetoku 

13.95% 

95.60 

99.79% 

85.59«S% 

Flat in Linear Space 

3.76% 

25.76 

14.88% 

88.81«»;|% 

Flat in Log Space 

8.54% 

58.51 

68.24% 

39.951^% 

Gaussian 

9.74% 

66.73 

79.08% 

79.34+20 ; 65% 

Specified Function 

9.05% 

62.00 

63.18% 

96.32™,% 
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Table 4: Summary of the results from simulations with A°g ak evolution. 


shown here are based on different E^ ak , as described in Sect. 


6.1 


The five samples 


rp src 
■^peak 

Distribution 

Detection 

Rate 

Prediction for 
Swift [yr” 1 ] 

KS-test significance 
for z distribution 

KS-test significance for 
peak flux distribution 

Modified Yonetoku 
Flat in Linear Space 
Flat in Log Space 
Gaussian 

Specified Function 

15.05% 

4.65% 

8.63% 

10.10% 

9.37% 

103.10 

31.86 

59.12 

69.19 

64.19 

98.62% 

31.39% 

59.15% 

76.84% 

57.44% 

69-34™% 

0 . 32 ^- 02 % 

47.881^% 

50-68™% 

89.00™% 


6.2. l?peak evolution 


As discussed in Sect. |3.2| , spectral evolution has been observed in many GRBs, and thus 
we implement an option to include A° aak evolution in our simulation. To see how the results 
might change if A°g ak evolution is included, we apply this option to the best-fit parameter 
set (Table |2j) with all five different Apeak distributions we test. 


Results with A° aak evolution included are summarized in Table |dj Similar to those 
samples presented in Table [3j results are from simulations of 10000 bursts and using the 
average number of 26884 active detectors from year 2005 to 2009. These simulations show 
that including A°g ak evolution can result in noticeable changes of the outcome, especially in 
the distributions of the burst characteristics. Although our best-fit sample (the one with 
a modified Yonetoku E^ ak distribution) remains good matches with both the redshift and 
peak-flux distributions, samples with other Ap a c ak distributions show that the resulting KS- 
test values can change considerably with A° aak evolution included. However, we noted that 
when actually plotting these distributions, the changes do not seem to be as significant as 
indicated by the KS-test significance levels in the tables. The general shapes and widths of 
the distributions remain similar with or without A° aak evolution. This is because the KS-test 
values are especially sensitive to the medium of the distribution. Therefore, a slight changes 
in the medium can lead to a major difference in the significance level. Moreover, if the 
uncertainties in the significance levels of the peak-flux distributions are taken into account, 
the data from real observations actually cannot distinguish the difference between the results 
with or without A° aak evolution. 


The sample using a flat E^ ak distribution in linear space (see description in Sect. 6.1) 
shows the most remarkable change in the KS-test significance level when A°j) ak evolution 
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is included. Therefore, Fig. 
evolution of this sample as an example of how the general shapes of the distributions change. 


14 


plots the peak-flux distributions with and without E °^ |k 


Peak Flux Distribution of the Mock-Triaaered Bursts 



Fig. 14. — Comparison of the peak-flux distributions with and without E °^ ® k evolution for 
the sample using a flat E^ k distribution in linear space. The result with E°^l k evolution is 
plotted as red bars. The distribution without E°^l k evolution is plotted as blue bars. Error 
bars along the y - axis show the statistical errors in each bin. Error bars along the x-axis 
represent the bin sizes. 


When including the E° ] ^ k evolution, the detection rates in these samples either stay 
similar or increase slightly compared to those without E°^l k evolution. A possible explanation 
is that as the burst spectra get softer, some of the bursts might have E°^l k fall into the BAT 
energy range, and thus become detectable, ffowever, due to the low number of samples with 
-^peak evolution in our simulations, we cannot exclude the possibility that this trend is a mere 
coincidence. 
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7. Selection biases of using single trigger criterion 

To explore possible selection biases when using a single trigger criterion, we use our 
best-fit sample (discussed in Sect. [5j) as an example and calculate the detection fraction for 
each criterion. The detection fraction here refers to the number of simulated bursts detected 
by each criterion divided by the total number of detections (when using all criteria). In 
other words, the detection fraction is the fraction of the mock-triggered bursts that would 
be detected if only using one criterion. Because most of the bursts are detected by more 
than one criterion, the detection fractions from all criteria do not add up to one. Results 
show that the detection fraction for each criterion varies from ~ 83% for the most efficient 
criterion (Trigger Criterion 7 ^ 334) to ~ 19% for the least efficient criterion (Trigger Criterion 
7 ^ 356). Therefore, adopting the complex trigger algorithm of BAT increases the detection 
rate by 1.2 to 5.3 times. The energy band 25 — 100 keV is the most sensitive energy range 
of BAT, which matches well with the peak of the BAT’s effective area. 



Fig. 15. — The peak-flux distributions of the simulated bursts triggered by the most efficient 
trigger criterion (Criterion 334) and the least efficient trigger criterion (Criterion 356). Error 
bars along the y - axis show the statistical errors in each bin. 
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Figure. [15] plots the peak-flux distributions of the simulated bursts triggered by the 
most efficient criterion (Trigger Criterion # 334) in red bars, and those triggered by the 
least efficient criterion (Trigger Criterion # 356) in blue bars. As one can see in the figure, 
the most efficient criterion is much more sensitive to the low-flux bursts than the least 
efficient criterion. Therefore, it is possible to miss some dim bursts when using only one 
single trigger criterion, and results in finding GRB characteristics that are biased toward 
the brighter end. The complex trigger algorithm of BAT generates more detections of the 
low-flux events, and thus improves the survey sensitivity. 


8. GRB detection fraction of the BAT-trigger algorithm 


To generally understand how the complex BAT-trigger algorithm affects burst detec- 
tions, we calculate the detection rate as a function of redshift. From redshift z = 0 to 10, 
we generate two hundred GRBs in each redshift bin with bin size of Az = 0.2. This is to 
make sure each redshift bin has enough bursts to create a statistically meaningful result. 
The redshift distribution of the bursts in each bin are uniformly assigned. Other than the 
redshift, all other burst characteristics (such as luminosity function, spectral distribution, 


etc), are the same as those found in our best-fit sample (Sect. 5.2) 


The results are plotted in Fig. 16. Panel (a) shows the fraction of detectable bursts as 


a function of redshift. Error bars along the y - axis show the statistical error in each bin (i.e. , 
\f~N, where N is the number of bursts in each bin). As expected, the detectable fraction is 
~ 1 at z ~ 0, and drops to approximately zero at high redshift. Above redshift z = 6, the 
detectable fraction is about 1 — 3%. Panel (b) presents the average flux of the detectable 
bursts in each redshift bin. The average flux of the detectable bursts decreases with redshift 
as anticipated, since bursts become dimmer when they are further away. At redshift z > 6, 
the average flux of the detectable bursts is around 10 -8 erg s -1 cm -2 , which is usually only 
detectable when the burst appears on-axis relative to the detector’s plane. There remain 
some statistical fluctuations in the plot, particularly at high redshift due to the small number 
of detections. 


9. Comparison with the star-formation rate 

As discussed in Sect, [lj long GRBs are related to at least a sub-class of core-collapse 
supernovae. Core-collapse supernovae are expected to directly trace the star-formation his- 
tory due to the short lifetimes of their massive progenitors. Thus, the cosmic long GRB rate 
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Fig. 16. — Panel (a): The fraction of detectable GRBs as a function of redshift. The result 
is based on an intrinsically flat GRB redshift distribution, with other GRB characteristics 


(e.g., luminosity and spectral distributions) from our best-fit sample (Sect. 5.2). Error bars 


along the y-axis show the statistical errors in each bin. Panel (b): The average flux of 
the detectable GRBs in the simulation, as a function of redshift. The adopted bin size is 
Az = 0.2. 


is expected to follow the shape of the cosmic SFR. Kistlcr et al. (2008a) noted that from the 


high-redshift GRB observations, there is an unexpected rise in the cosmic GRB rate at large 


redshift compared to that expected from previous SFR measurements. Yiiksel et al. (2008) 


use the high-redshift GRB measurements to calculate the high-redshift SFR, and conclude 


that SFR in the early universe might be larger than previously expected (Hopkins & Beacom 


2006). 
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The comparison between onr best-fit cosmic GRB rate and the shape of cosmic SFR 

Onr 


from both Hopkins & Beacom (2006) and Yiiksel et al. (2008) can be found in Fig. 17 


best-fit GRB rate is plotted as a red line. The green and blue lines in the figure show the 


GRB rates that strictly follow the shapes of the SFR from Hopkins & Beacom (2006) and 


Yiiksel et al. (2008), respectively. The normalizations of the green and blue lines come from 


fitting with real GRB observations by including luminosity evolution (see discussions in the 


following section, Sect. 10). 


The red shaded region in Fig. [17] shows the uncertainty of our best-fit GRB rate. We 
quantify the uncertainty by modifying the parameters in the GRB rate function (Eq. [Tj) 
around our best-fit set of parameters (Table [2]) until results no longer match well with 
observations. The shaded region in the figure indicates the parameter space that produces 
results which satisfy the following three criteria: (1) Matching with the observed peak- flux 
distribution with KS test significance level > 90%. (2) Matching with the observed redshift 
distribution with KS test signihcance level > 90%. (3) Prediction for the Swift detection 
rate within the range of 95 ± a/95 = 95 ± 10 events per year. Table [H] and [6] summarize 
the parameters and results for the lower and upper limit of the GRB rate shown in Fig. [17} 
Note that for the peak-flux distribution match, we require the KS test significance including 
the error bars exceeds 90%. Therefore, it is the upper limit of the KS test significance that 
needs to exceed 90%, as seen in the tables. 


Additionally, the constraining factor for this uncertainty region of the GRB rate turns 
out to be the prediction of the Swift detections (i.e. , the third criterion listed above), rather 
than the shapes of the functions. In other words, all the GRB rates within the red-shaded 
region produce decent matches with the shapes of the observed peak-flux and redshift dis- 
tributions. However, if we move the curve further away from the lower/upper limit of the 
red-shaded region, we will have the predicted detection rates lower/higher than the true rate 
of Swift , even though the shapes of the peak-flux and redshift distributions still match well 
with the observations. 



Horiuchi et al. 

(2009 

) quantify the uncertainty of the star-formation rate measurements 

summarized in 

Yiiksel et al. 

(2008 

) by taking into account the scatter of data. In order to 


compare with our best-fit sample and have a better idea of how well we can use the GRB rate 
to constraint the star-formation rate at high redshift, we plot this uncertainty as the blue 
shaded region in Fig. [17] with the same normalization factor used for the blue line. Note 
that the blue-shaded region is the uncertainty for the star-formation rate measurements, 
instead of the GRB rate. In other words, although the blue line can generate a good fit 
with the GRB observations using the current normalization if including luminosity evolution 
(see Sect. 10), there is no guarantee that GRB rates within the blue-shaded region can all 


41 


Table 5: Summary of the set of parameters for the lower limit of the GRB rate. Parameters 
for the luminosity function are the same as those shown in Table [2] 


-Rgrb(£ — 0) 
[Gpc -3 yr -1 ] 

Zl 

n\ 

n 2 

Detection 

Rate 

Prediction for 
Swift [yr -1 ] 

0.75 

3.60 

2.10 

-3.50 

18.70% 

83.85 


KS-test (D) for 
z distribution 

Significance for 
z distribution 

KS-test (D) for 
peak flux distribution 

Significance for 
peak flux distribution 

5.70 x 10 -2 

98.68% 

3.32+ 2 ;^ x 10 -2 

93.09+®^% 


Table 6: Summary of the set of parameters for the upper limit of the GRB rate. Parameters 
for the luminosity function are the same as those shown in Table [2] 


Rgrb(z — 0) 
[Gpc -3 yr -1 ] 

Zl 

n\ 

n 2 

Detection 

Rate 

Prediction for 
Swift [yr -1 ] 

1.01 

3.60 

1.95 

-0.00 

12.95% 

104.74 


KS-test (D) for 
z distribution 

Significance for 
z distribution 

KS-test (D) for 
peak flux distribution 

Significance for 
peak flux distribution 

5.41 x 10 -2 

99.43% 

5.12+ 3 ;]]® x 10 -2 

58.10^,372% 


generate good matches with observations. 


Results in Fig. 17 show a clear diversity in the shapes of the GRB rate versus redshift 
from the SFRs at z ~ 4. Both the GRB rate and the SFRs start decreasing beyond x ~ 4. 
However, the SFRs, even the one from Yiiksel et al. (2008), decline much faster than the 
GRB rate found in our best-fit sample. The SFRs from Hopkins & Beacom (2006) and 


Yiiksel et al. (2008) decrease at large redshift with a power-law index of ~ —8.0 and ~ —3.5, 


respectively, while the GRB rate found in this work decreases with a power-law index ~ —0.7. 


There are several possible explanations for this very high GRB rate at large redshift. 
First, this could suggest an even larger SFR at high redshift, which implies that most of the 
star formation activities at high redshift probably come from low- luminosity galaxies, and 
thus measurements of the SFR based on galaxy observations might underestimate the true 


Kistlcr et al. 

2009 

Jakobsson et al. 

2012; 

Tanvir et al. 

2012 

Kistler et al. 

2013 
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Trenti et al. 2013). Alternatively, a high GRB rate at early times can also be explained if the 
fraction of GRB-related supernovae changes as a function of redshift. That is, if there are 
more supernovae accompanied by GRBs at high redshift, one could get a higher GRB rate 


without adjusting the SFR. For example, Woosley & Heger (2012) suggest several collapsar 


models that can generate long gamma-ray transients, and state that these events are more 
frequent in the early universe. 

Another possibility would be the luminosity evolution. If we allow luminosity evolution 
in our simulation and generate more high-luminosity bursts at high redshift than at low 
redshift, this could create enough low-flux bursts without over-producing the detections at 
low redshift. In this case, we might not need a high GRB rate at large redshift to balance 
the low-flux ratio of those from the observed GRBs. Several studies have already suggested 


the possibility of redshift evolution of the GRB luminosity function (e.g., Salvaterra et al. 


2009, 

2012 

Virgili et al. 

2011 

Toma et al. 

2011 


will thus investigate this possibility in the following section. 


10. Possibility of Luminosity Evolution 


To explore whether the extreme excess of GRB rate at high redshift is not necessary 
when luminosity evolution is considered, we test several functions for luminosity evolution 
using the shapes of the previously reported SFRs as the intrinsic GRB rate. Our goal is to 
study whether it is possible to generate a sample that matches well with real observations, 
while having the GRB rate restricted to follow the current SFR measurements. Specifically, 


we perform the tests with two possible SFRs: (1) the commonly adopted SFR from Hopkins 
& Beacom (|2006), and (2) the one in Yiiksel et al.| (2008), which is similar to the first one 


but with high-redshift corrections using GRB detections. 


For simplicity, we made the characteristic luminosity L* in the luminosity function 
(Eq.§ change as a function of redshift z. Two different functional forms are tested in our 
simulations: (1) L* = A x z B , and (2) L* = A x logio(z). Again, all simulations with 
luminosity evolution are based on 26884 active detectors, which is the average number from 
year 2005 to year 2009. 


Results show that the second form (L* oc logi 0 (z)) can create mock-triggered samples 
that match well with the observations, with some adjustments of the parameters in the 
luminosity function. Table [7] and [8] summarize the GRB characteristics that generate a good 


match with observations when adopting luminosity evolution and the SFR from [Hopkins 


& Beacom (2006) and Yiiksel et al. (2008), respectively. The KS-test values for these two 
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Fig. 17. — Comparison between the cosmic GRB rate from our best-fit sample (red line; 


Sect. 5.2) and the cosmic GRB rates that follow strictly the shapes of the SFRs from Hopkins 


& Beacom (2006) (green line) and Yiiksel et al. (2008) (blue line). The GRB rates that trace 
the SFRs can generate results that match well with observational data if luminosity evolution 


is included (see discussion in Sect. 10). 


samples comparing with the observed redshift and peak-flux distributions are also given in 
the tables. As expected, a more severe luminosity evolution is needed if assuming the SFR 


from Hopkins & Beacom (2006), because this rate decreases more rapidly at high redshift 


than that from Yiiksel et al. (2008). 


The green and blue lines in Fig. 17 plot these two GRB rates with characteristics 
summarized in Table [7] and |8j respectively. As discussed before, the GRB rate shown as the 


green line strictly follows the shape of the SFR from Hopkins & Beacom (2006), while the 
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Table 7: Summary of the parameters and results of the best-fit sample with luminosity 

evolution using the SFR from Hopkins & Beacom (2006). 


Rgrb(z = 0)[Gpc 3 yr 1 

Functional Form 

1.07 

SFR in Hopkins & Beacom (|2006j) 


u 

[erg s _1 ] 

A 

X 

y 

zpsrc 

■^peak 

distribution 

Detection 

Rate 

Prediction for 
Swift [yr -1 ] 

1Q51.00 

2.0 

-0.20 

-2.00 

Flat in Log Space 

22.28% 

95.27 


KS-test (D) for 
z distribution 

Significance for 
z distribution 

KS-test (D) for 
peak flux distribution 

Significance for 
peak flux distribution 

5.96 x 10~ 2 

97.26% 

2.42+ 3 ; 21 x 10" 2 

98.76^8% 


Table 8: Summary of the parameters and results of the best-fit sample with luminosity 

evolution using the SFR from Yiiksel et ah (2008). 


Rgrb(z = 0)[Gpc 3 yr 1 

Functional Form 

1.11 

SFR in 

Yiiksel et al. 

(2008 

) 


u 

[erg s _1 ] 

A 

X 

y 

zpsrc 

■^peak 

distribution 

Detection 

Rate 

Prediction for 
Swift [yr -1 ] 

1Q51.00 

1.9 

-0.20 

-2.00 

Flat in Log Space 

21.01% 

95.11 


KS-test (D) for 
z distribution 

Significance for 
z distribution 

KS-test (D) for 
peak flux distribution 

Significance for 
peak flux distribution 

6.13 x IQ’ 2 

96.48% 

2.34 + 3 ;™ x 10- 2 

99.19^3%% 


GRB rate in blue traces the shape of the SFR from Yiiksel et al. (2008). Our search suggests 
that these two GRB rates can only match well with real GRB observations if luminosity 
evolution is included. 


45 


11. Discussions and Conclusions 


We developed a program that simulates the complex trigger algorithm adopted by BAT. 
We used this program to search for a cosmic GRB rate and luminosity function that generates 
a mock-triggered sample with characteristics that match well with observations. Our results 
suggest that the BAT’s complex trigger algorithm increases the detection rate by ~ 1.2 to 5.3 
times higher than that using a single flux threshold. Therefore, adopting the complex trigger 
algorithm of BAT improves the chance of triggering bursts with low fluxes. As a result, we 
need an intrinsic GRB sample that is on average dimmer than previously expected, in order 
to avoid over-producing the number of detections and to match with real Swift observations. 
This can be achieved by either adding more bursts with lower luminosities, increasing the 
number of bursts at high redshift, or both. 

According to all the parameter sets we tried, generating more bursts with lower lu- 
minosities in the intrinsic GRB luminosity function has a side effect of creating too many 
detections at low redshift, and thus resulting in a distribution that does not match well with 
the redshift distribution of the real BAT-detected GRBs. Therefore, adding more bursts at 
large redshift provides a way to increase the number of dim bursts without over-producing 
low-redshift observations, and thus generates a mock-triggered sample that matches well 
with both the redshift and the peak-flux distributions of real observations. 


By adopting the complex trigger algorithm of BAT, the best- fit sample we found suggests 
the possibility of an intrinsic GRB rate that contains even more bursts at higher redshift than 
previous expectations (Yiiksel et al. 2008 Butler et al. 2010 Wanderman & Piran 2010). 
This result implies that either (1) the SFR in the early universe is even higher than that 
inferred by previous GRB studies, and hence a majority of star formation might happen in 
low- luminosity galaxies, or (2) some redshift evolution effects, such as luminosity evolution 
or an evolution in the GRB-to-supernovae ratio, need to be taken into account. 


Theoretical studies suggest that the GRB characteristics are likely to be different in the 


early universe (e.g., Salvaterra et al. 2009 Virgili et al. 2011; Toma et al. 2011; Woosley & 


Heger 2012). Therefore, at least some redshift evolution in the GRB luminosity functions are 


expected. We thus examine the possibility of including luminosity evolution to reduce the 
intrinsic GRB rate at high redshift and maintain a good match with observations. We found 
that if we assume the shape of the SFR from Yiiksel et al. (2008), the characteristic luminosity 
L* in Eq. [3] needs to evolve as L* = 1.9 x logi 0 (z) in order to generate a mock-triggered sample 
that matches well with both the observed redshift and peak-flux distributions. If we assume 
the shape of the SFR from Hopkins & Beacom (2006), which predicts an even lower GRB 
rate at high redshift, a more severe luminosity evolution with L* = 2.0 x log 10 (z) is needed 
to produce a good match with observations. 
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Besides the GRB luminosity, other GRB characteristics might also evolve with redshift. 
Theoretical studies suggest that low-metallicity population III stars produce GRBs that have 


much longer durations than regular ones 

(Fryer et al. 

2001 

Komissarov & Barkov 

2010 

Meszaros & Rees 

2010 

Toma et al. 

2011 

Wooslcy & Heger 

2012 

). For example, 

Gendre 


et al. (2013) propose that the ultra-long GRB 111209A, which had a prompt emission that 
lasted around 1.5 x 10 4 seconds, resulted from a low-metallicity blue supergiant star and 
resembled more the population III star explosions. Therefore, one might expect GRBs in 
the early universe to have longer durations. The BAT trigger algorithm, particularly the 
rate trigger, is relatively insensitive to these ultra-long bursts due to their slow-changing 
light curve. Our trigger simulator shows that an ultra- long burst like GRB 111209A (with 
similar pulse duration, E iso , and spectral parameters) can be detected by image trigger out 
to redshift z ~ 1.17 if the event happens on-axis relative to the detector plane (Grid ID 
= 17), and z ~ 0.3 if the burst is far off-axis (Grid ID = 14). The rate trigger criteria 
can only detect such bursts within z ~ 0.17 even if the burst appears on axis. Note that 
in this particular simulation, we still assume a flat background level throughout the burst 
duration, which is likely not to be true when the event lasts for a few hours. Therefore, such 
long-duration GRBs are even harder to be detected by the BAT in practice, and thus might 
introduce further uncertainty of the intrinsic GRB rate at high redshift. 

The best-fit sample from our simulations suggests the GRB rate in the local universe 
to be ~ 0.84 Gpc~ 3 yr _1 , which is in general agreement (within a factor of 2) with other 
observations and studies (e.g., Schmidt 2001 Guetta & Della Valle 2007 Liang et al. 2007 


Pelangeon et al. 2008} Wanderman & Piran 2010). Note that this is the rate of GRBs 


beamed toward us, as we do not consider bursts that are pointed away from us, due to the 
large uncertainty in the beaming factor (see discussion in Sect [Tj) . In addition, the best-fit 
sample suggests that BAT should detect ~ 96 GRBs per year, which is consistent with the 
real BAT detection of ~ 95 bursts per year averaged from 2005 to 2009 (Sakamoto et al. 
2011). According to our best-fit sample, ~ 1% of all detections (i.e. , ~ 1 burst per year) 


come from redshift z > 6. 

Moreover, the GRB rate from our best result gives a total number of 4568^ 3 42 9 GRBs 
per year that are beamed toward us in the whole universe. The errors here are calculated 
using the lower and upper limits of the GRB rate shown in Fig. [XT] (the red-shaded region). 
To have a better idea how this number is compares to the cosmic star-formation history 
and the core-collapse supernova rate, we perform a simple order-of-magnitude calculation, 


as described below. Using the current star-formation rate measurements (e.g., Hopkins & 


Beacom i 

CO 

o 

o 

Yiiksel et al. 2008 

) and the commonly-adopted modified Salpeter Initial Mass 

Function 

(Salpeter 

1955 

Baldry & Glazebrook 2003 

), one can estimate the total core-collapse 


supernova rate in the whole universe to be ~ 10 8 per year (for example, see calculation in 
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Lien & Fields 2009 

. Observations suggest that ~ 25% of all the core-collapse supernovae 

are Type Ibc event 

Li et al. 2011), and only ~ 1% of all Type Ibc supernovae are related to 

GRBs (e.g., 

Berger et al. 

2003 

). An or der-of- magnitude calculation then gives ~ 2.5 x 10 5 


GRB per year in the whole universe. If we assume the beaming factor to be ~ 50 (i.e., 1/50 
GRBs are pointed at us; Guetta et al. 2005), there will be ~ 5000 GRBs per year that are 
beamed toward us, which is consistent with our result. 


Furthermore, for all the different -E^k distributions we tried, it seems to be difficult to 
generate results that match well with the observed distribution without assuming some 
kind of intrinsic correlation between 77 p g C ak and the burst energy output, such as the L pea k . 
Similar conclusions have been drawn by many groups that study possible correlations in GRB 
characteristics (e.g., Ghirlanda et al. 2012 Shahmoradi 2013). Better U// ak measurements 


will be crucial to verify this conclusion. 


We presented here a GRB rate and luminosity function that generates a mock-trigger 
sample that matches well with observation. To confirm whether there exist equally good or 
better fits other than the one we presented, and also to quantify the probability of whether 
the good match is due to a real physical solution or just a coincidental match, a complete 
Monte Carlo search would be required. Therefore, a significant decrease of the simulation 
time will be important to further pin down the cosmic GRB rate, improve our understanding 
of the GRB characteristics, and quantify their uncertainties. 


Observationally, increasing the number of GRB detections with redshift measurement is 
essential in reducing the uncertainties in the observed distributions and better constraining 
the GRB properties. Moreover, better knowledge of the biases in GRB follow-up observations 
and redshift determination is critical to understand the completeness of the observational 
GRB sample. The rapidly growing number of GRB detections by Swift and the effort in 
measuring GRB characteristics through multi-wavelength observations will certainly improve 
our understanding of GRB physics, and will make them even better probes of stellar evolution 
and star-formation history out to the early universe. 
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